#############################################
# File to compute and plot substantive      #
# results for:                              #
# "The Arc of Modernization: Economic       #
# Structure, Materialism, and the           #
# Onset of Civil Conflict"                  #
#                                           #
# J. Tyson Chatagnier and Emanuele Castelli #
#                                           #
# Forthcoming in Political Science Research #
# and Methods                               #
#                                           #
# File created 31 March 2016                #
#############################################

rm(list = ls())

library(foreign)
library(ggplot2)

##############################
# Write function for 3D plot #
##############################

VarPlot3d <- function(x.vals,
                      y.vals,
                      p.vec,
                      cols.vec    = c("darkblue",
                                      "lightblue"
                                      ), 
                      par.box     = NULL,
                      perspective = TRUE,
                      screen.x    = -60,
                      screen.y    = 0,
                      screen.z    = -150,
                      xlab        = "X Variable",
                      ylab        = "Y Variable",
                      zlab        = "Probability \n of War",
                      ...
                      ) {
                        # This creates a 3d plot, varying both x and y.
                        # x.vals is our first variable
                        # y.vals is our second variable
                        # p.vec is our vector of probabilities (pre-calculated)
                        # cols.vec is an optional vector of color arguments
                        # for the wireframe plot; must be of length 2
                        # The remaining arguments are optional arguments to 
                        # the wireframe function
  require(lattice)
  vars.df <- data.frame(x = x.vals,
                        y = y.vals,
                        z = p.vec
                        )
  newcols <- colorRampPalette(c(cols.vec[1],
                                cols.vec[2]
                                )
                               )
  wireframe(z ~ x * y,  
            data        = vars.df, 
            drape       = T, 
            screen      = list(z = screen.z, 
                               y = screen.y,
                               x = screen.x
                               ), 
            col.regions = newcols(length(p.vec)), 
            par.box     = par.box,
            perspective = perspective,
            xlab        = xlab, 
            ylab        = ylab, 
            zlab        = zlab,
            ...
            )
}

################
# Read in data #
################

cw <- read.dta("cw_replication_data.dta")

####################################
# Get coefficients and covariance  #
# matrices from the Stata analysis #
####################################

b <- read.table("sppo_coefs_final.txt")
vcov <- read.table("sppo_vcov_final.txt")

coef.names <- names(b)

b <- as.numeric(b)

names(b) <- coef.names

vcov <- as.matrix(vcov)

##############################
# Make profile for predicted #
# probabilities.             #
##############################

mod.seq <- seq(from       = min(cw$mod,
                                na.rm = TRUE
                                ),
               to         = max(cw$mat,
                                na.rm = TRUE
                                ),
               length.out = 40
               )

mat.seq <- seq(from       = min(cw$mat,
                                na.rm = TRUE
                                ),
               to         = max(cw$mat,
                                na.rm = TRUE
                                ),
               length.out = 40
               )

mod.mat <- expand.grid(mod.seq,
                       mat.seq
                       )

pred.df <- data.frame(mod       = mod.mat[, 1],
                      mat       = mod.mat[, 2],
                      oil       = 0,
                      gdp       = mean(cw$gdp,
                                       na.rm = TRUE
                                       ),
                      growth    = mean(cw$growth,
                                       na.rm = TRUE
                                       ),
                      polity    = -1,
                      constant  = 1,
                      pop       = mean(cw$pop,
                                       na.rm = TRUE
                                       ),
                      mtn       = 0,
                      ef        = mean(cw$ef,
                                       na.rm = TRUE
                                       ),
                      parity    = mean(cw$parity,
                                       na.rm = TRUE
                                       ),
                      sixties   = 0,
                      seventies = 0,
                      eighties  = 0
                      )

pred.df$p2 <- pred.df$polity ^ 2

pred.cons <- rep(1,
                 times = nrow(pred.df)
                 )


r.sq <- with(pred.df, 
             cbind(mod,
                   mat,
                   oil,
                   sixties,
                   seventies,
                   eighties,
                   gdp,
                   growth,
                   polity,
                   p2,
                   constant
                   )
             )

r.w <- with(pred.df, 
            cbind(pop,
                  mtn,
                  ef,
                  parity
                  )
            )

r.acq <- matrix(pred.df$constant)

g.w <- with(pred.df, 
            cbind(mod,
                  oil,
                  gdp,
                  p2,
                  pop,
                  sixties,
                  seventies,
                  eighties,
                  constant
                  )
            )

b.rsq <- b[1 : 11]
b.rw <- b[12 : 15]
b.racq <- b[16]
b.gw <- b[17 : 25]

pb <- pnorm(g.w %*% b.gw)
pa <- pnorm(pb * (r.w %*% b.rw) + (1 - pb) * (r.acq %*% b.racq) - r.sq %*% b.rsq)

pw <- pa * pb

#################
# Plot Figure 4 #
#################

pdf("mod_mat_3d.pdf")
VarPlot3d(x.vals   = mod.mat[, 1],
          y.vals   = mod.mat[, 2],
          p.vec    = pw,
          xlab     = "Modernity",
          ylab     = "Materialism",
          cols.vec = c("blue", 
                       "aliceblue"
                       ),
          screen.x = 320, 
          screen.y = 0, 
          screen.z = 160,
          main = "Effects of Modernity and Materialism on Civil Conflict"
          )
dev.off()

##########################
# Results for real cases #
##########################

#####################
# Results for Chile #
#####################

chl.df <- na.omit(cw[cw$ccode == 155, ])

r.sq <- with(chl.df, 
             cbind(mod,
                   mat,
                   oil,
                   sixties,
                   seventies,
                   eighties,
                   gdp,
                   growth,
                   polity,
                   p2,
                   1
                   )
             )

r.w <- with(chl.df, 
            cbind(pop,
                  mtn,
                  ef,
                  parity
                  )
            )

r.acq <- matrix(1,
                nrow = nrow(r.sq)
                )

g.w <- with(chl.df, 
            cbind(mod,
                  oil,
                  gdp,
                  p2,
                  pop,
                  sixties,
                  seventies,
                  eighties,
                  1
                  )
            )

b.rsq <- b[1 : 11]
b.rw <- b[12 : 15]
b.racq <- b[16]
b.gw <- b[17 : 25]

pb <- pnorm(g.w %*% b.gw)
pa <- pnorm(pb * (r.w %*% b.rw) + (1 - pb) * (r.acq %*% b.racq) - r.sq %*% b.rsq)

pw <- pa * pb

chl.p <- data.frame(year = chl.df$year,
                    prob = pw
                    )

act.conf <- chl.p[chl.df$onset == 1, ]

###############################
# Plot left panel of Figure 5 #
###############################

pdf("chile_sppo.pdf")
  ggplot(chl.p,
         aes(x = year,
             y = prob
             )
         ) + 
         geom_point(size = 2.5) +
         geom_smooth(colour = "navyblue",
                     size = 1.1
                     ) +
         geom_point(data = act.conf,
                    aes(x = year,
                        y = prob),
                    shape  = 8,
                    colour = "red",
                    size   = 2.5
                    ) +
         xlab("Year") +
         ylab("Probability of Civil War Onset") +
         theme_bw()
dev.off()

######################
# Results for Mexico #
######################

mex.df <- na.omit(cw[cw$ccode == 70, ])

r.sq <- with(mex.df, 
             cbind(mod,
                   mat,
                   oil,
                   sixties,
                   seventies,
                   eighties,
                   gdp,
                   growth,
                   polity,
                   p2,
                   1
                   )
             )

r.w <- with(mex.df, 
            cbind(pop,
                  mtn,
                  ef,
                  parity
                  )
            )

r.acq <- matrix(1,
                nrow = nrow(r.sq)
                )

g.w <- with(mex.df, 
            cbind(mod,
                  oil,
                  gdp,
                  p2,
                  pop,
                  sixties,
                  seventies,
                  eighties,
                  1
                  )
            )

b.rsq <- b[1 : 11]
b.rw <- b[12 : 15]
b.racq <- b[16]
b.gw <- b[17 : 25]

pb <- pnorm(g.w %*% b.gw)
pa <- pnorm(pb * (r.w %*% b.rw) + (1 - pb) * (r.acq %*% b.racq) - r.sq %*% b.rsq)

pw <- pa * pb

mex.p <- data.frame(year = mex.df$year,
                    prob = pw
                    )

act.conf <- mex.p[mex.df$onset == 1, ]

###############################
# Plot left panel of Figure 6 #
###############################

pdf("mexico_sppo.pdf")
  ggplot(mex.p,
         aes(x = year,
             y = prob
             )
         ) + 
         geom_point(size = 2.5) +
         geom_smooth(colour = "navyblue",
                     size = 1.1
                     ) +
         geom_point(data = act.conf,
                    aes(x = year,
                        y = prob),
                    shape  = 8,
                    colour = "red",
                    size   = 2.5
                    ) +
         xlab("Year") +
         ylab("Probability of Civil War Onset") +
         theme_bw()
dev.off()

